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NOMENCLATURE 
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T;  T ; T 
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constant  in  burning  rate  law,  r = apn 

Eo  - w-  ln'd) 

specific  heat  of  solid  and  gas,  respectively 
Eg/R0^So»  [n"d] ; activation  energy  of  surface  reaction 
defined  in  Eq.  (lc),  [n-d] 

Q /c  (T  - T ) , [n-d] 

thermal  conductivity  of  solid  and  gas,  respectively 
defined  in  Eqs.  (5)  and  (6),  respectively 
exponent  in  Eq.  (9) 

mass  flow  rate  = pgr  [see  Eq,  (2a)] 
pressure  exponent  in  burning  rate  law,  r = apn 
pressure  exponent  in  Eq.  (13) 
pressure;  p/p0  [n-d] 
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heat  release  per  unit  mass  for  surface  reaction  and 
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surface  regression  rate;  r/rQ,  [n-d] 
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temperature;  initial  surface  temperature;  cold-solid 
temperature 

distance  coordinate,  normal  to  propellant  surface 
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I . INTRODUCTION 


Theoretical  efforts  to  predict  gun  interior  ballistics  are  faced 
with  the  formidable  problem  of  analyzing  a three-dimensional,  viscous, 
heat  conducting,  chemically  reacting,  two-phase  flow  field.  Because 
of  these  difficulties,  recent  work  has  concentrated  on  a simpler,  one- 
dimensional, inviscid  model  which  then  requires  correlation  formulas 
to  determine  drag  losses  and  heat  transfer  rates  between  the  gas  and 
solid  particles.  Equally  important  is  an  'a  priori'  description  of  the 
behavior  of  the  solid  propellant  burning,  rate.  All  current  gun  interior 
ballistics  computer  programs  (1-6)  rely  on,  and  are  sensitive  to,  the 
steady  state  burning  rate  "law",  r = ap  , possibly  modified  with  an  ero- 
sive burning  contribution.  However,  this  relationship  is  assumed  to  be 
valid  in  a highly  transient  pressure  environment.  If  the  burning  rate 
is  not  properly  represented,  a successful  gun-code  prediction  of  maxi- 
mum pressure  and  muzzle  velocity  may  indicate  that  other  adjustable 
parameters  in  the  interior  ballistics  theory  have  artificially  compen- 
sated. 


Gough , P.  S.  and  Z wants,  F.  J. , "Theoretical  Model  for  Ignition  of 
Gun  Propellant 3"  Report  SRC-R-67,  Space  Research  Corp.  (December  1972); 
also  Goughj  P.  S.3  " Fundamental  Investigation  of  the  Interior  Ballis- 
tics of  Guns 3"  Report  SRC-R-74,  Space  Research  Corp . (August  1974). 

2 

Frier,  H. , Van  Tassell3  W. , Ragan , S.,  and  VerShaw,  J.  T.  3 "Model  of 
Gun  Propellant  Flame  Spreading  and  Combustion ,"  BRL-CR-147  March  3 
AD  #9188421;  also  Frier 3 H. s " Predictions  of  Flamespreading  and 
Pressure  Wave  Propagation  in  Propellant  Beds3"  AAE  7S-63  University 
of  Illinois  (July  1975). 

3 

Kuo3  K. , Vichnevetsky 3 R. 3 and  Summer field3  M.  , "Theory  of  Flame  Front 
Propagation  in  Porous  Propellant  Charges  Under  Confinement 3 " AIAA 
Journal 3 Vol  11 3 No.  Z3  pp  444-451  (April  1973). 

4 

Fisher 3 E.  B.  3 and  Trippe3  A.  P. s "Development  of  a Basis  for  Accep- 
tance of  Continuously  Produced  Propellant3  " Report  VQ-5S163-D-lt 
Calspan  Corp.  (November  19-73). 

5 

East3  J.  L. 3 and  McClure3  D.  R.3  "Projectile  Motion  Predicted  by  a 
Solid/Gas  Flow  Interior  Ballistic  Model3  " 10th  JANNAF  Combustion 
Meeting3  CPIA  publication  243  (August  1973). 

ft 

Baer3  P.  G. 3 and  Frankie 3 J.  M. , "The  Simulation  of  Interior  Ballistic 
Performance  of  Guns  by  Digital  Computer  Program 3"  BRL  Report  1183 
(December  2982),  AD  #299980. 
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Experimental  verification  of  the  instantaneous  propellant  burning 
rate  during  a rapid  pressure  change  is  a difficult  task  and,  at  the 
present  time,  is  unavailable.  Similarly,  the  time-dependent  numerical 
solution  of  a comprehensive  solid  propellant  combustion  model  which  in- 
cludes solid  and  gas-phase  finite-rate  chemical  reactions  has  yet  to  be 
accomplished.  An  alternative  is  to  examine  the  thermal-wave  combustion 
theories  (based  on  the  quasi-steady  flame  assumption)  which  were  origi- 
nally developed  for  stp<iy-*of  combustion- instability  and/or  extin- 
guishment in  solid  propellant  rocket  engines.  It  is  the  purpose  of  this 
report  to  investigate  the  behavior  of  three  such  combustion  models  in  a 
rapidly  increasing  pressure  field  typical  of  a large-caliber  gun.  The 
three  models  selected  are: 

1.  KTSS  (Krier,  Tien,  Sirignano,  Summerfield)  - Ref,  7. 

2.  Levine/Culick  - Ref.  8. 

3.  Kooker/Zinn  - Ref.  9. 

The  primary  objective  is  to  obtain  the  numerical  solution  for  instan- 
taneous propellant  burning  rate  for  a series  of  example  cases  where  the 
imposed  pressure  variation  at  the  edge  of  the  flame  zone  is  prescribed 
by  an  experimental  pressure-time  history  from  the  105mm  Gun  M68,  These 
results  can  then  be  compared  to  the  steady  state  expression,  r = apn.  A 
secondary  objective  is  to  evaluate  the  capabilities  of  four  different 
numerical  integration  schemes  which  are  used  to  generate  the  time-depen- 
dent solution  to  these  solid  propellant  combustion  models. 

II.  ANALYSIS 

A.  Derivation  of  Combustion  Models 

The  derivation  of  each  of  the  three  models  is  based  on  a common  set 
of  simplifying  assumptions.  It  is  assumed  that  (1)  the  solid  propellant 
combustion  process  can  be  described  in  a single  spatial  dimension,  (2) 
the  unburned  propellant  is  homogeneous,  incompressible,  and  inert,  (3) 
the  solid  is  converted  into  gas  by  a global  pyrolysis  reaction  which 
occurs  at  an  infinitesimally  thin  interface,  and  (4)  the  gas-phase  flame 


Krier,  H. , T'ien,  J.  S.,  Sirignano , W.  A.,  and  Simnerfield,  M.,  " Non- 
steady Burning  Phenomena  of  Solid  Propellants:  Theory  and  Experiments," 

AIAA  Journal , Vol.  6,  No . 2,  pp  278-28S  (February  1968). 

8 

Levine,  J.  N .,  and  Culiak,  F . E.  C .,  "Nonlinear  Analysis  of  Solid 
Rocket  Combustion  Instability,"  AFRPL-TR-7 4-45  (October  1974). 

9 

Kooker,  D.  E.,  and  Z inn,  B.  T. , "Numerical  Investigation  of  Nonlinear 
Axial  Instabilities  in  Solid  Rocket  Motors BRL-CR-141  (March  1974), 

AD  0776954,  see  also  AIAA  Paper  73-1298. 
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is  quasi-steady  and  remains  anchored  to  the  interface.  As  a direct  re- 
sult of  these  assumptions,  all  three  models,  when  written  for  the  coordi- 
nate system  shown  in  Fig.  1,  reduce  to  the  following  initial/boundary 
value  problem  describing  the  thermal-wave  in  the  unburned  solid  propellant 


Figure  1 - Combustion  Model  Configuration 


0 

T 

+ R0-0  =0  (-  « < n < 0) 

n nn 

(1) 

initial 

condition 

0(n,  0)  = eRn 

(a) 

boundary 

( 0(-«  , x)  = 0 

(b) 

conditions 

( en(0,  x)  = G (0S,  R,  P) 

(c) 

with  R = R(0g) 

(d) 

where  the  particular  form  of  G and  R depends  on  the  model  in  question. 
The  quasi-steady  flame  assumption  is  incompatible  with  a proper  des- 
cription of  the  ignition  event.  Hence,  the  initial  condition  in  the 
present  application  is  taken  to  be  the  thermal-wave  corresponding  to  an 
assumed  steady-state  burning  rate.  Equation  (lb)  of  the  boundary  con- 
ditions enforces  the  cold-solid  temperature  at  infinity,  while  Eq.  (lc) 
is  the  result  of  an  energy  balance  at  the  interface  and  includes  the 
contribution  of  the  gas-phase  flame.  The  point  of  departure  of  the 
three  models  is  the  flame  analysis  and  the  postulated  global  pyrolysis 
reaction  [Eq.  (Id)]. 

When  the  low-velocity  flame  zone  is  assumed  to  be  quasi-steady  and 
one-dimensional,  the  equations  of  global  continuity  and  momentum  yield 
the  almost  trivial  results. 


11 


pu  = p r = m = ih(t] 
g g s 


(a) 


(2) 


p = p(t).  (b) 

If  the  coupling  with  the  species  continuity  equations  is  essentially 
ignored,  the  energy  equation  can  be  written  in  the  form. 


. dT  d ,,  dT\ 

mcp  dx  * dx  Ckg  dx3  = V (3) 

where  the  right-hand  side  is  the  product  of  the  heat  released  in  the 
flame  reaction (s)  and  the  mass  production  rate  of  product  species.  The 
gradient  boundary  condition  [Eq.  (lc)]  is  a statement  of  the  energy 
balance  across  the  interface,  x = 0,  which  can  be  written  symbolically 
as 

/ heat  flux  into  /heat  released  in  \ /heat  flux  from  ( 

\ solid  propellant/  (surface  reaction/  jgas-phase  flame/ 

The  "heat  feedback"  from  the  flame  zone  follows  from  the  integral  of 
Eq.  (3),  evaluated  at  the  interface.  Before  the  integration  can  be 
carried  out,  the  spatial  distribution  of  must  be  specified  'a  priori'. 

If  the  problem  were  posed  correctly,  the  spatial  distribution  in 
question  could  easily  be  computed  after  obtaining  the  simultaneous 
solution  of  the  energy  equation  and  the  species  continuity  equations. 

To  assume  the  distribution  of  Q^co  is  known  beforehand  implies,  in  effect, 
that  the  solution  has  already  been  determined.  Without  resolving  this 
contradiction,  all  three  models  employ  the  distribution. 


Qfm 


a constant 

0 


0 < x < xf 
x > x^ 


(5) 


where  Xf  is  the  location  of  the  edge  of  the  flame  zone.  Another  possi- 
bility would  be 


QfW  = K2  5(x  - xf)  (6) 

where  <$  is  the  Dirac  delta  function.  Obviously,  in  both  cases,  the 
"constant"  K will  change  with  time.  It  is  sometimes  claimed  that  in- 
tegrating the  energy  equation  under  the  constraint  of  Eq.  (5)  means 
that  the  flame  analysis  approximates  a diffusion-controlled  process. 
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whereas  use  of  Eq.  (6)  (flame-sheet)  approximates  a chemical  kinetic 
rate-controlled  process.  On  this  basis  alone,  the  analysis  which 
follows  from  Eq.  (S)  is  said  to  apply  to  a composite  propellant  while 
the  analysis  following  from  Eq.  (6)  should  apply  to  a double-base  pro- 
pellant. Considering  the  complexity  of  the  actual  solid  propellant 
combustion  process  and  the  degree  of  simplicity  and  approximation  in- 
volved in  the  one-dimensional  models,  it  is  felt  that  such  a distinction 
is  highly  speculative.  With  due  caution,  the  present  analysis  will 
assume  the  three  models  based  on  Eq.  (5)  represent  the  combustion  of  a 
triple-base  propellant  (M30)  and  interpret  the  results  as  possible 
trends . 

The  different  form  of  the  gradient  boundary  condition  [Eq.  (lc)] 
in  each  model  can  be  traced  to  the  method  of  evaluating  Kj  in  Eq.  (5) 
after  obtaining  the  integral  of  Eq.  (3).  In  the  KTSS  model,  when  the 
characteristic  combustion  time  is  assumed  to  be  a function  of  pressure 
alone,  the  integral  of  the  energy  equation  leads  to  the  conclusion  that, 

qg  = ®(p)/r  (7) 

where  qg  is  the  heat  feedback  from  the  gas-phase  flame  to  the  interface. 
$(p)  is  determined  by  assuming  its  functional  dependence  on  pressure  in 
the  unsteady  environment  is  the  same  as  during  steady  state  combustion, 
i.e.,  <t(p)  = r qg.  The  final  result  for  the  gradient  boundary  condition 
is  (see  Ref.  7 for  details) 

0 (0,t)  = RH  + [P2n  (Pn/m  - H)]/R  (8) 

where  R = 0 m.  (9) 

s 

The  power  law  dependence  of  R on  surface  temperature  is  used  as  an 
approximation  to  an  Arrhenius  reaction  rate.  The  appearance  of  the 
pressure  exponent  "n"  in  Eq.  (8)  is  the  result  of  assuming  that  the 
steady  state  burning  rate  of  the  propellant  is  r = apn.  Thus,  if  the 
pressure  is  held  constant,  the  time-dependent  model  will  predict  a 
value  of  burning  rate  identical  to  this  steady  state  law  when  all  tran- 
sient portions  of  the  solution  have  died  away.  Of  course,  for  a given 
value  of  pressure  in  an  arbitrary  time-dependent  situation,  the  heat 
feedback  from  the  gas-phase  flame  will  not  necessarily  equal  its  steady 
state  value  at  that  pressure. 

In  the  Levine/Culick  model,  the  value  of  Kj  (and  its  dependence  on 
pressure)  is  also  determined  without  specifying  the  actual  reaction 
scheme  in  the  flame  zone.  Combining  the  steady-state  temperature 
gradient  in  the  solid  at  x = 0,  the  integral  of  Eq.  (3)  evaluated  at 
x = 0,  and  the  steady-state  energy  balance  at  the  interface  yields  an 
expression  for  w,  viz. 
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By  assuming  that 


ih  = psapn  = Bgp  s exp  (-ES/R0T$)  (11). 

allows  m and  Ts_to  be  eliminated  from  Eq.  (10),  leaving  only  the  pressure 
dependence  of  u>.  Then  in  the  spirit  of  the  quasi-steady  flame  assump- 
tion, the  equation  governing  u>  is  found  by  replacing  the  steady  state 
value  of  pressure  with  its  instantaneous  value.  After  considerable 
algebra  (see  Ref.  8) , the  final  result  for  the  gradient  boundary  con- 
dition is 


0n(o,T)  = R[H  + (0s  - 1)  Cl  - ^)]  + P2n  j (1  - H)  + 


!e  E r . 1 

c A L n - n 

i + _L_ — lnP 


1]  i /R 


(12) 


ns  A(es  - i) 

where  R = P exp  [ — 


1 + f (0S  - « 


(13) 


Although  somewhat  disguised,  Eq.  (13)  is  the  unsteady  equivalent  to  the 
Arrhenius  expression  in  Eq.  (11);  for  all  computations  shown  in  the 
present  paper,  ng  = 0.  Similar  to  Ref.  7,  the  Levine/Culick  model  will 
also  predict  a value  of  burning  rate  identical  to  ap  1 when  the  pressure 
is  held  constant  and  the  transient  portion  of  the  solution  has  vanished. 


In  the  Kooker/Zinn  model,  the  flame  region  is  assumed  to  be  con- 
trolled by  the  simple  one-step  global  reaction. 

Oxidizer  + Fuel  Products  (14) 


which  implies 

8 

0)  ~ p • 


(15) 


Although  Eq.  (14)  is  a global  reaction  as  opposed  to  an  elementary  one, 
it  is  not  improbable  that  8,  the  "reaction  order",  could  be  2 as  used 
in  Ref.  9.  The  present  computations  also  employed  the  value  1.7.  Under 
the  assumption  that  the  reaction  occurs  at  a known  flame  temperature 


and  that  the  reaction  order  is  specified,  the  integral  of  Eq.  (3)  leads 
to  the  result, 


Cd 

0n(O,T)  = R[H  + (0  - 1)  Cl  - ^)]  + Z (16) 

s 

with  R given  by  Eq.  (13)  when  ns  = 0.  The  unknown  multiplicative  con- 
stant associated  with  Eq.  (15)  is  determined  by  the  interface  energy 
balance  at  the  initial  condition,  i.e., 

Z = 1 - H (17) 

where  the  reference  state  in  the  nondimens ionalizat ion  is  taken  to  be 
the  given  initial  conditions.  Since  a particular  burning  rate  law  has 
not  been  used  in  the  derivation  of  Eq.  (16),  it  is  of  interest  to  com- 
pute the  pressure  dependence  of  the  steady  burning  rate  inherent  in 
the  model.  The  example  in  Fig.  2 (based  on  the  propellant  parameters 
of  Table  1)  shows  the  results  of  constant  pressure  calculations  for^ 
several  values  of  pressure  (hence  the  symbols  in  Fig,  2),  assuming  B = 
2,  The  plots  of  lnRs  vs.  InP  are  not  straight  lines,  but  gentle 
curves  concave  upward.  For  lower  values  of  surface  heat  release  (H  < 
0.7  in  this  example^  the  curvature  is  slight;  thus  the  model  closely 
approximates  r = ap  . For  larger  values  of  H,  the  curvature  is  more 
pronounced;  hence  the  model  predicts  a substantial  increase  in  burning 
rate  "exponent"  with  increasing  pressure.  In  this  example,  both  com- 
putations (H  = 0.70  and  0.90)  were  unrealistically  forced  through  the 
same  initial  conditions.  This  is  responsible  for  the  prediction  that 
the  larger  value  of  H leads  to  a lower  overall  burning  rate.  However, 
the  important  point  is  that  a very  simple  flame  model  can  predict  an 
increase  in  burning  rate  exponent  with  increasing  pressure  as  observed 
experimentally  for  double  and  triple-base  propellants. 

B.  Numerical  Solution  Methods 


Predictions  of  the  time-dependent  burning  rate  for  each  combustion 
model  require  the  solution  to  Eq.  (1)  using  the  appropriate  interface 
boundary  condition  [Eq.  (8)  for  KTSS;  Eq.  (12)  for  Levine/Culick;  Eq. 
(16)  for  Kooker/Zinn] , An  exact  analytical  solution  is  not  possible 
since  the  primary  equation  and  the  boundary  conditions  are  coupled  and 
nonlinear.  Because  of  the  potential  difficulty  introduced  by  the  non- 
linearity, the  present  study  examined  four  different  numerical  solution 
techniques,  each  of  which  has  been  used  in  previous  investigations  to 
solve  combustion  problems  of  this  type.  The  objective  is  to  see  if 
the  desired  solution  is  independent  of  the  method  used  to  obtain  it. 


Figure  2.  Steady  State  Burning  Rate  Predicted  by  Kooker/Zinn  Model 


Each  numerical  technique  is  expressed  in  the  following  notation. 

The  superscripts  n and  n+1  refer  to  the  current  time  level  (where  the 
solution  is  known)  and  the  new  time  level,  respectively.  The  subscripts 
j-1,  3 » 3+1  refer  to  an  arbitrary  interior  mesh  point,  j,  and  its  two 
adjacent  neighbors.  For  reasons  to  be  discussed  below,  the  mesh  point 
distribution  is  non-uniform  and  hence  An j = nj  - nj-i  and  An j + 1 = 
r\.+j  - r\.  are  not  necessarily  equal.  For  convenience,  define 


D0 . 
J 


2 2 2 2 
An.  ©.  i + (An  . , - An  .)  0.  - An  . , 0.  , 

...J J.+I __J*1 _n  1L  ^ ill 


An . An . , (An . + An . , ) 
3 3+1  3 3+1 


(18) 


= 0 at  point  j , 


DD0  = 

3 An.An.+1(An.  + Anj+1) 


(19) 


50  at  point  j . 
nn 

Then  the  four  numerical  techniques  when  applied  to  Eq.  (1)  can  be 
written  as, 

2 12 
(A)  Explicit  [0(At,  An  ),  stable  for  Ax<  ^ An 

r 1 v ’ max'  2 mnJ 


0.n+1  - 0.n 

-2 T i-  + Rn  (D0.n)  - DD0,n  = 0 . 

At  3 j 


(20) 


(B)  Full  Implicit  [0(Ax,  An  ),  stable  for  any  Ax) 

max 

e.n+1  - 0.n 

-J L.  + Rn+1  (D0.n+1)  - DD0.n+1  = 0 . 

At  k 3 ' 3 


(21) 


2 2 

(C)  Crank-Nicolson  Implicit  [0(At  , An  ),  stable  for  any  At] 

max 


e.n+1  - © n 


At 


j-  + J j Rn+1  (D0jn+1)  - DD0jn+1 


(22) 


+ Rn(D0^n)  - DDOj11  J = 0 


(D)  Invariant  Imbedding  Implicit  [0(At,  An  ) stable  for  any  At] 

max 

See  Appendix. 
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The  coupling  between  the  unknowns  at  j-1,  j,  and  j+1  in  Methods  (B)  and 
(C)  leads  to  a tridiagonal  matrix  which  is  inverted  with  the  standard 
Thomas  Algorithm. 

All  the  methods  generate  a finite  difference  equation  (FDE)  which 
is  only  an  approximation  to  the  true  partial  differential  equation  (PDE). 
The  relationship  between  the  two  is  based  on  Taylor  series  expansions 
and,  using  Method  (B)  as  an  example,  can  be  written  as, 


At 


FDE  = PDE  - — ^ ®TT+  [RAynj+1  + 2(An.  - Anj+1ncenm/6) 


+ [RAHj  Ar)j+j  (Anj  + 1 - Arij)  - 2{br\*  - An^An.^  + An^j)] 


(0  /24)  . 

nnn  ' 


C23) 


For  equally- spaced  mesh  points,  this  reduces  to  the  more  familiar  form. 


FDE  = PDE  - 0 + (2R0  - 0 )(An2/12)  . (24) 

2 tt  v nnn  nnn  ' 

2 ■ 

Based  on  Eq.  (24),  the  formal  order  of  accuracy  is  0(At)  and  0(An  ). 
However,  unless  the  derivative  terms  are  uniformly  small,  the  formal 
argument  can  be  misleading.  Since  the  present  combustion  problem  in- 
volves a thermal  profile  which  may  have  very  steep  gradients  near  the 
interface,  the  spatial  mesh  must  be  closely  spaced  in  this  region  (i.e.. 

An  <<  1).  Otherwise  the  truncation  error  will  swamp  the  desired  solution. 
Although  no  single  mesh  system  will  be  optimal  for  all  problems.  Method 
(D)  was  used  to  determine  a 231  grid  point  system  [21  points,  -0.01  < 
n < 0;  30  points,  -0.10  < n < -0.01;  50  points,  -2.00  < n < -0.10;  130 
points,  -15,00  < n < -2.00]  such  that  the  numerical  solution  for  an 
extreme  combustion  problem  was  independent  of  further  rearrangement  and 
additional  points. 


The  problem  of  nonlinearity  must  also  be  treated  with  care.  Each 
numerical  method  effectively  linearizes  Eq.  (1)  in  some  manner.  Further- 
more the  interface  boundary  condition,  which  must  be  applied  at  the  new 
time  level  n+1,  is  a nonlinear  function  of  the  surface  temperature  at 
n+1,  This  requires  an  iteration  procedure  where  the  integration  of 
Eq.  (1)  and  the  application  of  the  boundary  condition  are  updated 
cyclically.  Convergence  is  assumed  when  the  difference  between  the  R 
used  in  Eq.  (1)  and  the  R which  follows  from  the  boundary  condition  is 
less  than  a specified  tolerance  (in  this  study,  normally  10"^  using 
double  precision  computation).  A crucial  point  concerning  the  size  of 
the  time  step  At  enters  here.  Although  Methods  (B)  through  (D)  are 
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theoretically  stable  for  an  unrestricted  time  step,  the  iteration  pro- 
cedure is  not.  For  those  combustion  problems  where  the  surface  tempera- 
ture is  changing  rapidly  with  time,  the  time  step  must  be  appropriately 
restricted  to  keep  the  numerical  solution  from  diverging.  No  mathe- 
matical proof  is  offered  to  support  this  conclusion.  However,  in  the 
present  investigation,  the  time  step  was  continuously  monitored  and 
appropriately  reduced  if  the  number  of  cycles  required  for  convergence 
of  the  iteration  exceeded  five. 

C.  Model  Equation 

The  limiting  accuracy  of  the  numerical  solution  computed  by  Methods 
(A)  through  (D)  can  be  estimated  with  an  example  denoted  as  the  model 
equation.  An  exact  analytical  solution  is  available  (see  p.  389,  Ref. 
10)  since  both  the  equation  and  the  boundary  conditions  are  linear. 

Using  the  present  notation,  this  system  can  be  written  as 


0 

T 

initial 

+ R0  - 0 =0  [t  > x , -®  < n < 0] 

n nn  L o'  - J 

(25) 

condition 

G(n,  tq)  = 0o(n) 

(a) 

boundary 

0(-«,t)  * 0 

(b) 

conditions 

T 

0n(O,x)  = H[0(O,  x)  + - ^T-] 

(c) 

so  “ 


where  R and  H are  constants. 

The  example  plotted  in  Fig.  3,  with  R=4.0,  H=5.0,  and  xo=0.01  is 
viewed  as  a reasonable  test  of  the  methods  even  though  it  simulates  a 
very  large  surface  heating  rate.  The  predicted  numerical  solution, 
based  on  the  same  grid  mesh  system  used  in  the  combustion  problems,  and 
the  exact  solution  are  indistinguishable  on  the  scale  of  Fig.  3.  The 
maximum  error  occurs  at  the  surface  n=0,  and  was  within  1%  for  Method  (D) . 
Furthermore,  using  Method  (D) , the  magnitude  of  the  error  at  a given 
time  level  decreased  as  the  size  of  the  time  step  At  decreased.  Method 
(A)  may  be  equally  accurate  for  this  linear  problem  but  the  diffusion 
stability  restriction  on  the  time  step  is  so  small  that  the  method  is 
not  competitive  economically.  For  this  reason,  it  was  eliminated.  For 
Methods  (B)  and  (C) , the  maximum  surface  temperature  error  was  within 
2%,  but  the  expected  convergence  as  the  time  step  was  decreased  did  not 
always  materialize.  This  can  be  traced  to  errors  in  the  tridiagonal 
matrix  inversion  used  by  both  Methods  (B)  and  (C) . Using  a uniform  mesh 
system  as  an  example,  the  diagonal  element  of  the  matrix  is  given  by 


10 

Car  slaw,  H.  S. , and  Jaeger , J.  C. , Conduction  of  Heat  in  Solids , 
Oxford  Press , 2nd  edition  (1959). 
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(26) 


1 + 


2Ax 

An2 


while  the  off-diagonal  elements  are 


RAx  + At 

2An  . 2 
An 


(27) 


For  time  steps  at  or  below  the  diffusion  stability  limit  and  with 
R <<  2/An,  diagonal  dominance  is  assured.  However,  when  using  large 
time  steps  and  nonlinear  spacing  (variable  An) , the  absolute  value  of 
the  diagonal  term  can  change  by  orders  of  magnitude  from  one  corner  of 
the  matrix  to  the  other.  This  can  create  the  same  type  of  errors  in 
the  matrix  inversion  process  as  does  loss  of  diagonal  dominance.  These 
errors  can  easily  overpower  the  temporal  truncation  error  in  the  finite 
difference  equation,  leading  to  the  result  that  a decrease  in  Ax  does 
not  produce  the  expected  increase  in  accuracy.  Fot  the  example  shown 
in  Fig.  3,  however,  the  loss  in  accuracy  is  minimal. 

On  the  basis  of  the  linear  model  equation,  it  is  concluded  that 
Methods  (B)  through  (D)  yield  nearly  equal  results  with  acceptable 
accuracy.  The  results  from  the  linear  model  equation  are  not  intended 
to  be  conclusive  proof  of  accuracy  for  the  combustion  models  to  follow, 
but  rather  an  estimate  of  the  expected  lower  error  bound.  It  should 
be  stressed  that  the  influence  of  the  nonlinearity  in  the  three  com- 
bustion models  has  not  been  assessed  in  this  exercise. 

III.  RESULTS 

A.  Burning  Rate  Response  to  a Sudden  Compression 

The  purpose  of  this  investigation  is  to  compute  the  burning  rate 
response  to  a rapidly  increasing  pressure  field  which  simulates  the  gun 
combustion  chamber  environment.  A similar  problem  is  posed  in  Ref.  7 
where  one  area  of  concern  is  the  behavior  of  a rocket  engine  during  a 
sudden  compression.  Several  numerical  computations  in  Ref.  7 are  based 
on  the  pressure  profile  shown  in  Fig.  4,  which  is  given  by 

P(x)  = 1.0  + AP  [1.0  - exp(-Bx)]  . (28) 

When  AP  = 2.5  and  8 = 1.0,  this,  expression  approximates  an  experimental 
test  run  reported  in  Ref.  7.  Equation  (28)  is  not  a good  representation 
of  a typical  gun  pressure  time-history.  However,  it  is  relevant  to  the 
gun  problem  if  viewed  as  a large  amplitude  pressure  wave  suddenly 
passing  over  a propellant  grain  which  is  burning  at  steady  state.  Such 
a situation  could  conceivably  develop  during  a malfunction.  For  this 
reason,  any  results  which  follow  from  Eq.  (28)  are  of  interest  to  the 
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present  analysis.  A set  of  computations  showing  the  sensitivity  of  the 
KTSS  model  to  changes  in  the  value  of  surface  heat  release  (H)  is  re- 
produced here  as  Fig.  S (see  Fig.  9 in  Ref.  7).  The  unbounded  burning 
rate  predicted  for  H = 0.80  is  one  of  several  calculations  used  as 
evidence  of  "intrinsic  burning  rate  instability"  or  "runaway".  A 
similar  runaway  condition  was  reported  in  Ref.  8 as  the  result  of  a 
sinusoidal  pressure  oscillation. 

Because  of  the  adverse  influence  that  burning  rate  runaway  would 
have  in  a gun  system,  the  decision  was  made  to  examine  each  model  in 
detail  to  determine  the  cause  of  the  unbounded  result  and  when  it  might 
occur.  A first  step  was  to  use  Method  (D)  to  recompute  the  KTSS  re- 
sults for  H = 0.80  shown  in  Fig.  5.  The  recomputed  results  are  shown 
in  Fig.  6;  no  runaway  condition  is  indicated  and  the  burning  rate  re- 
turns to  the  same  value  as  H = 0.75  at  a nondimensional  time  of  approxi- 
mately 1.2.  Further  computations  of  this  case,  however,  showed  that 
Method  (D)  will  also  predict  a divergent  solution  if  the  time  step 
restriction  imposed  by  the  iteration  procedure  is  removed  or  a suffi- 
cient number  of  mesh  points  near  the  propellant  interface  are  eliminated. 

To  see  if  "runaway"  exists  under  more  extreme  conditions,  the 
value  of  surface  heat  release  (H)  was  increased  to  0.90  and  the  response 
to  the  same  exponential  pressure  profile  was  computed.  The  results  (Method 
D)  shown  in  Fig.  7 contain  sharp  finite-amplitude  "spikes",  similar  to 
those  reported  in  Ref.  9.  An  important  difference,  however,  is  that  the 
burning  rate  response  in  Fig.  7 is  due  to  a monotonically  increasing 
pressure  field  and  not  a sinusoidal  pressure  oscillation.  The  behavior 
of  the  thermal  wave  during  a spike  is  the  same  as  that  discussed  on  page 
24  pertaining  to  the  results  in  Fig.  14.  Although  these  results  appear 
unconventional  and  may  suggest  that  the  combustion  model  is  oversimplified, 
an  extensive  rechecking  of  the  numerical  computation  indicates  that  Fig. 

7 is  the  actual  solution  to  the  equations.  To  determine  if  the  finite- 
amplitude  spikes  are  model  dependent,  the  Kooker/Zinn  combustion  model 
(with  H = 0.9  and  8 = 1.7)  was  driven  by  the  same  monotonically  in- 
creasing pressure  field  given  by  Eq.  (28).  The  numerical  solution 
(Method  D)  shown  in  Fig.  8,  exhibits  the  same  qualitative  behavior  as 
the  KTSS  model.  A strong  resemblance  to  this  type  of  burning  rate 
phenomenon  has  been  reported  in  the  Russian  literature  (Ref.  11)  for 
a related  problem.  Reference  11  considers  "gas-less"  combustion  of  a 
solid  propellant  (or  explosive)  which  burns  by  means  of  a single,  in- 
depth,  irreversible,  condensed-phase  reaction.  The  equations  solved 
are  identical  to  those  used  by  Bradley  (Ref.  12,  with  a - 1)  in  his 
study  of  radiant  ignition  of  a reactive  solid.  Since  the  model  does 


11 

Shkadinskii,  K.  G, , Khaikin , B.  I. , and  Merzhanov3  A.  G. , "Propagation 
of  a Pulsating  Exothermic  Reaction  Front  in  the  Condensed  Phase 3 " 
Fizika  Goreniya  i Varyva3  Vol  ?3  No.  1,  pp  19-28  (1971). 
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Bradley , H.  H. , Jr. , "Theory  of  Ignition  of  a Reactive  Solid  by  Con- 
stant Energy  Flux3"  Comb  Sci  & Tech>  Vol  2,  pp  11-20  (1970). 
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Figure  5.  Burning  Rate  Response  of  KTSS  Model  (Original  Results) 
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Figure  6.  Burning  Rate  Response  of  KTSS  Model  (^©computed  Results) 


Figure  7.  Burning  Rate  Response  of  KTSS  Model 
With  H=0, 9 
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not  consider  regression  of  a gas/solid  interface,  the  burning  rate  quoted 
in  Ref.  11  is  the  velocity  of  the  reaction  front  maximum  as  it  moves  in- 
to the  unreacted  solid.  For  certain  values  of  heat  release,  the  velocity 
of  the  reaction  front  pulsates  (see  Fig.  5 in  Ref.  11)  in  a manner  simi- 
lar to  the  present  results  for  surface  regression  shown  in  Figs.  7 and  8. 

It  is  also  of  interest  to  compare  the  predictions  of  all  three  com- 
bustion models  under  identical  input  conditions  when  driven  by  the  same 
exponential  pressure  profile  [Eq.  (28)].  Figure  9 shows  the  results 
(Method  D)  of  such  a comparison  when  H = 0.8,  n = 0.67,  and  the  remaining 
parameters  are  those  of  Table  I. 


Table  I.  Assumed  Propellant  Parameters 


3 

pg  = 1.54  gm/cm 

cs  = 1.55  J/gm-°K  [0.37  cal/gm-°K] 

kg  = 0.0031  J/cm-sec-°K  [7.3  x 10-4  cal/cm-sec-°K] 

-3  2 

=>a  = 1.28  x 10  cm  /sec 


= 298°K 


= 0.0628  MJ/mole  [15  kcal/mole] 

= 1.72  J/gm-°K  [0.41  cal/gm-°K] 

= 0.00156  J/cm-sec-°K  [3.72  x 10 
= 0.75  cm/sec  \ 

=6.9  MPa  [1000  psia]  [ 


(m  = 6.2,  KTSS) 

4 cal/cm-sec-°K] 


steady  state  initial  conditions 


= 623°K 


Qualitatively,  the  burning  rate  response  of  each  model  is  the  same. 
The  large  value  of  peak  burning  rate  predicted  by  the  KTSS  model  can 
be  traced  to  the  assumption  that  Cp  = cs,  which  removes  the  damping 
term,  R(0S  - 1) (1  - Cp/cs) , from  the  gradient  boundary  condition. 
Since  it  could  easily  be  replaced  in  the  model,  this  is  considered  a 
minor  point  in  the  present  investigation. 
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Figure  9.  Comparison  of  Model  Predictions  for  Burning  Rate  Response 

to  KTSS  Pressure  Profile 


Two  important  observations  can  be  made  from  the  results  displayed 
in  Figs,  6-9,  First,  the  interval  of  time  equal  to  a/r  2 (At=1,  non- 
dimensional  ly)  , often  quoted  in  theoretical  analyses  as  the  characteris- 
tic response  time  of  the  thermal  wave  in  the  unburned  solid,  is  of 
limited  value  when  estimating  the  time  for  a significant  change  in  the 
burning  rate.  Secondly,  the  instantaneous  propellant  burning  rate 
cannot  be  characterized  as  a function  of  dlnP/dt,  To  dramatize  the 
latter  point,  consider  the  following  expression  for  instantaneous 
burning  rate  derived  recently  by  Krier  (Ref.  13)  from  a perturbation 
analysis  of  the  KTSS  Model.  In  the  notation  of  the  present  report, 

Rto  = Rs[l  + C29) 

R 


where  Y (P) 


a - - pn/m) 
v ro 

[Pn/m(2-^)  - 2H) 


and  P = P(x)  . 

Using  Rs  = Pn  and  P(t)  given  by  the  exponential  profile  in  Eq.  (28),  the 
numerical  predictions  of  Eq.  (29)  are  compared  in  Fig.  10  to  the  full 
numerical  solution  (Method  D)  of  the  KTSS  combustion  model  for  the  same 
set  of  parameters.  Certainly  for  this  example,  n<3  correlation  exists 
between  the  two  results  in  the  time  interval  0 - 1.0. 


B.  Burning  Rate  Response  in  a Simulated  Gun  Combustion  Chamber 

A typical  gun  propellant  is  assumed  to  be  represented  by  the 
parameters  given  in  Table  1.  Also  required  is  a value  for  the  surface 
heat  release  (H)  associated  with  the  global  pyrolysis  reaction.  This 
is  possibly  the  most  difficult  parameter  to  estimate  since  a single  re- 
action which  converts  solid  to  gas  does  not  exist  in  the  actual  com- 
bustion process.  The  present  investigation  will  rely  on  experimental 
studies  such  as  Ref.  14  in  which  estimates  are  made  of  the  amount  of 
heat  released  in  the  surface  zone  of  a double-base  propellant,  with 
and  without  catalysts,  under  steady  state  conditions.  The  results  of 
Ref.  14  are  reproduced  here  as  Fig.  11  which  indicates  that  a value 
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Figure  10.  Comparison  of  Equation  29  with  the  Numerical  Solution 

(Method  0)  of  the  KTSS  Model 


of  approximately  418  J/gm  (100  cal/gm)  would  be  compatible  with  the 
assumed  initial  conditions  of  Table  1.  This  implies  a range  of  values 
for  H between  0.8  and  0.9.  Other  investigators  (Refs.  15-17)  have 
inferred  values  of  H from  0.5  to  0.85  for  double-base  propellants. 
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Figure  11.  Surface  Heat  Release  vs  Burning  Rate  (Ref.  14) 

The  pressure  variation  to  be  imposed  at  the  edge  of  the  flame 
zone  in  the  combustion  models  is  obtained  from  an  experimental  firing 
of  the  105mm  Tank  Gun  M68  and  is  illustrated  in  Fig.  12.  The  record 
from  this  mid-chamber  gauge  indicates  an  ignition  and  flame-spreading 
delay  of  about  14  ms.  For  the  purpose  of  the  numerical  solution,  a 
time  origin  (t=0)  is  established  at  the  point  where  the  gun  pressure 
is  equal  to  6.9MPa  (1000  psi) . From  this  origin,  approximately  1.3 
nondimensional  time  units  (~3ms)  are  required  to  reach  the  peak 
pressure  of  393  MPa  (57,000  psi). 

An  example  of  the  burning  rate  response  to  the  above-mentioned 
gun  pressure-time  history  is  shown  in  Figs.  13  and  14  for  the  Levine/ 
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Figure  12.  105mm  Gun  Pressure  Time-History 
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Figure  13.  Influence  of  a Change  in  Surface  Heat  Release  on  Burning  Rate  Response  to  Gun 

Pressure  Time-History 
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Figure  14.  Burning  Rate  Response  to  Gun  Pressure  Time-History 


Culick  model  using  Method  D.  The  results  are  plotted  as  R/Rs,  where 
Rs  is  the  steady  state  value  of  burning  rate  based  on  the  instantaneous 
value  of  pressure  and  the  pressure  exponent,  n = 0.67.  Note  that,  at 
first,  R/Rs  falls  below  one,  indicating  that  the  thermal  wave  which 
controls  burning  rate  does  not  respond  instantaneously  to  the  increasing 
pressure.  The  figures  clearly  show  the  sensitivity  of  the  combustion 
model  as  H increases.  For  a value  of  H = 0.9,  the  numerical  solution 
indicates  a maximum  dynamic  overshoot  of  approximately  20.  To  see  if 
this  behavior  follows  from  unrealistic  variations  of  the  temperature 
distribution  in  the  solid  propellant,  the  thermal  profile  for  three 
points  in  time  (as  denoted  by  the  symbols  in  Fig.  14)  is  plotted  in 
Fig.  15.  The  curve  for  t = 0 is  the  initial  thermal  wave  at  steady 
state  which  decays  to  1/e  of  the  surface  value  in  a distance  of 
ct/rQ  (or  n = 1) . When  the  burning  rate  is  maximum,  the  distribution 
shows  a surface  temperature  of  880°K  and  a very  steep  gradient  in  the 
region  close  to  the  interface.  This  is  the  result  of  a heating  rate 
at  n = 0+  (surface  reaction  plus  flame  reaction)  which  has  been  greater 
than  the  rate  at  which  heat  can  be  conducted  into  the  unburned  solid 
without  raising  the  surface  temperature.  At  the  maximum  burning  rate, 
this  imbalance  reverses.  The  burning  rate  then  declines  as  a result 
of  the  combined  effects  of  (1)  rapid  heat  conduction  into  the  unbumed 
solid  due  to  the  steep  thermal  profile  near  the  surface,  and  (2)  the 
large  regression  rate  which  is  continuously  moving  lower  temperature 
propellant  closer  to  the  surface.  In  this  example  as  well  as  those 
shown  in  Figs.  7 and  8,  the  declining  burning  rate  is  accelerated  by 
the  exponential  dependence  of  the  surface  heat  release  term  on  pro- 
pellant surface  temperature  (i. e. , the  RH  term).  The  profile  at  t = 
0.106  shows  that  the  maximum  thermal  gradient  is  located  in  the  in- 
terior region.  At  this  lower  value  of  burning  rate,  the  gradient  at 
the  surface  is  much  less  than  its  previous  maximum  at  x = 0.080,  and 
hence  the  surface  heat  conduction  rate  is  much  lower.  However,  the 
rate  of  heat  release  in  the  flame  zone  continues  to  increase  with  the  in 
creasing  pressure.  Thus,  a new  dynamic  cycle  begins.  The  multiple 
spike  pattern  in  Figs.  7 and  8 is  created  in  the  same  way. 

The  influence  of  the  numerical  method  of  integration  (Methods  B,  C, 
D)  on  the  final  burning  rate  prediction  is  examined  in  Fig.  16  using  the 
Levine/Culick  combustion  model  with  H = 0.80  and  n = 0.67  as  an  example. 
The  results  show  a close  agreement  except  near  the  peak  burning  rate. 

In  the  vicinity  of  the  peak.  Method  D predicts  a larger  value  than 
Methods  (B)  and  (C)  which  are  coincident.  After  considerable  study,  it 
was  determined  that  the  difference  is  attributable  to  errors  in  the 
matrix  inversion  process  used  by  both  Methods  (B)  and  (C) . For  these 
implicit  solutions  where  the  time  step  is  much  larger  than  the  explicit 
diffusion  stability  limit,  a loss  of  diagonal  dominance  can  occur  when 
the  regression  rate  R is  large  (see  Eqs.  26  and  27).  In  a special  ex- 
ercise, the  time  step  At  and  the  grid  spacing  An.  near  the  propellant 
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Figure  15.  Temperature  Profiles  ip  Unburned  Solid  Propellant 
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Figure  16.  Comparison  of  Burning  Rate  Response  Predicted  by  Three  Numerical 

Integration  Methods 


interface  were  manipulated  to  maintain  diagonal  dominance  as  the  re- 
gression rate  became  large;  the  numerical  result  for  peak  burning  rate 
then  matched  the  Method  (D)  result  to  within  a fraction  of  a percent. 

This  verifies  the  hypothesis  but  the  required  manipulation  is  uneconomi- 
cal. It  should  be  emphasized  that,  at  all  times,  the  numerical  solutions 
obtained  with  Methods  (B)  and  (C)  were  stable  and  smooth.  Prof.  Gino 
Moretti's  well-known  statement  about  the  dangers  of  "confusing  smooth- 
ness with  accuracy"  is  definitely  applicable  here.  For  reasons  beyond 
the  scope  of  this  report.  Method  (D)  (invariant  imbedding)  is  not  sub- 
ject to  the  diagonal  dominance  limitation  and  is  considered  the  superior 
method,  in  both  accuracy  and  economy. 

The  influence  of  a change  in  the  steady  state  pressure  exponent  n 
is  demonstrated  in  Fig,  17  with  the  Levine/Culick  model,  again  using 
Method  D.  The  larger  value  of  n leads  to  approximately  the  same  rela- 
tive burning  rate  response,  but  at  an  earlier  time.  The  response  of 
all  three  combustion  models  to  the  gun  pressure-time  history  is  shown 
in  Fig.  18  using  the  larger  value  of  n.  Again  the  predictions  are 
qualitatively  the  same;  in  all  the  models,  the  dynamic  burning  rate  is 
ultimately  controlled  by  the  response  of  the  thermal  wave  in  the  un- 
bumed  solid. 

C.  Burning  Rate  Response  in  a Simulated  Closed  Bomb 

A simple  problem  was  devised  to  estimate  the  effect  of  dynamic 
burning  rate  on  closed-bomb  measurements.  The  propellant  is  assumed 
to  be  burning  at  steady  state  until  t = 0 when  it  is  suddenly  confined. 
The  imposed  pressure  is  then  given  by 

P(t)  = 1.0  + /qT  Rdt  . (30) 

The  numerical  predictions  (Method  D)  following  from  the  Levine/Culick 
model  (with  the  same  parameter  values  as  in  Fig.  18)  and  the  steady 
state  "model"  R = Pn  are  compared  in  Fig.  19.  A dynamic  overshoot  of 
approximately  3 occurs  at  low  pressure  (early  time)  but  the  visible 
effect  disappears  rapidly  thereafter.  If  the  results  below  35  MPa 
(5  kpsi)  were  discarded  when  analyzing  this  simulated  closed-bomb 
experiment,  the  dynamic  effects  would  go  undetected. 

CONCLUSIONS 

1.  "Intrinsic  burning  rate  instability"  (or  runaway)  is  a numerical 
difficulty  and  is  not  a solution  to  the  thermal-wave  combustion  models. 
The  nonlinear  models  admit  large-amplitude,  sharp  "spikes"  in  burning 
rate  as  a response  to  a monotonically  increasing  pressure  field,  but 
the  numerical  solutions  remain  finite. 
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Figure  17.  Influence  of  a Change  in  Pressure  Index  "n” 
on  Burning  Rate  Response  to  Gun  Pressure  Time-History 
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Figure  18.  Comparison  of  Model  Predictions  for  Burning  Rate 
Response  to  Gun  Pressure  Time-History 
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Figure  19.  Simulated  Closed-Bomb  Results  Using  the  Levine/Culick  Model  with  H=0.8,  n=0.8 


2.  When  the  combustion  models  are  driven  by  a gun  pressure-time 
history  or  a simulated  closed-bomb  environment,  the  combustion  response 
at  low  pressure  (early  time)  can  be  much  greater  than  the  "steady  state" 
prediction  provided  by  apn,  while  the  response  at  high  pressure  returns 
asymptotically  to  apn.  For  the  closed-bomb  problem  considered  here, 
the  dynamic  response  of  the  propellant  would  go  undetected  if  the 
pressure  record  below  35  MPa  (5  kpsi)  were  discarded. 

3.  Given  the  major  assumption  of  a quasi -steady  flame  region,  all 
three  combustion  models  predict  nearly  the  same  burning  rate  response 
regardless  of  the  differences  in  the  analysis  of  the  gas-phase  flame 
zone.  For  these  models,  the  dominant  influence  is  the  response  of  the 
thermal  wave  in  the  unburned  solid. 

2 

4.  The  interval  of  time  equal  to  a/r0  , often  quoted  in  theoretical 
analyses  as  the  characteristic  response  time  of  the  solid  propellant, 
can  be  misleading  when  used  to  estimate  the  time  for  a significant 
change  in  the  burning  rate.  The  present  results  show  that  substantial 
changes  in  burning  rate  can  occur  in  a time  interval  which  is  an  order 
of  magnitude  smaller  than  the  characteristic  response  time. 

5.  Simplified  expressions  for  instantaneous  propellant  burning  rate, 
derived  from  a given  combustion  model  using  linearized  analysis  (such 
as  Ref.  13),  can  yield  erroneous  trends  when  compared  to  the  numerical 
solution  of  the  complete  model. 

6.  Based  on  the  three  combustion  models  investigated  here,  the  dynamic 
burning  rate  is  very  sensitive  to  the  amount  of  heat  released  in  the 
surface  reaction.  Future  modeling  efforts  should  concentrate  on  a 
detailed  description  of  the  physical  and  chemical  changes  which  occur 
at  or  near  the  propellant  surface. 

Important  Implication  of  Conclusions 

The  present  investigation  has  demonstrated  that  three  thermal- 
wave  solid  propellant  combustion  models,  based  on  the  assumption  of  a 
quasi-steady  flame  zone,  predict  the  possibility  of  rapid,  large- 
amplitude  dynamic  burning  rate  excursions  from  the  steady  state  law, 
apn.  Except  for  the  simulated  closed-bomb  example,  these  calculations 
have  assumed  that  the  pressure  variations  and  the  propellant  burning 
rate  are  uncoupled,  i.e.  the  pressure  field  is  simply  imposed  at  the 
edge  of  the  flame  zone.  Of  course,  in  the  actual  gun  combustion 
chamber,  these  two  phenomena  are  coupled.  In  close  analogy  to  the 
combustion  instability  problem  in  a rocket  engine,  the  pressure 
disturbance  and  the  dynamic  burning  rate  response  it  creates  could  be 
"in-phase"  with  each  other,  where  one  will  reinforce  the  other.  The 
study  of  rocket  instability  in  Ref.  9 found  that  under  in-phase  con- 
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ditions  a dynamic  burning  rate  response  of  3-4  would  sustain  a longi- 
tudinal pressure  disturbance  in  the  rocket  engine;  it  should  be  noted 
here  that  the  continuous  outflow  through  the  .nozzle  represents  a sub- 
stantial loss  mechanism  to  the  disturbance  flow  field  in  the  chamber. 
If  the  sum  of  the  loss  mechanisms  inherent  in  the  gun  combustion 
chamber  flow  field  are  the  same  order  of  magnitude  as  a rocket  nozzle, 
then  a similar  behavior  would  be  expected.  Hence,  once  a pressure 
disturbance  is  created  in  a gun  system,  the  dynamic  burning  rate  could 
easily  be  responsible  for  driving  and  sustaining  a large  amplitude 
pressure  wave  in  the  combustion  chamber. 
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APPENDIX 


Application  of  Invariant  Imbedding 

Following  Ref.  Al,  Invariant  Imbedding  [Method  (D)]  is  applied  to 
a second-order  linear  partial  differential  equation  in  the  following 
manner.  Consider  the  system  (for  an  arbitrary  dependent  variable,  u) 

u.  + Ru  - u =0  [-x  < x < 0]  (Al) 

tXXX  1 O " “ J v ' 

initial  u(x'°>  * uoW  ta) 

condition 

Iu(-xQ,t)  = 0 (b) 

■ 

ux(0,t)  = G (u, t) . (c) 

Using  the  method-of-lines,  let 

ut  s (u  - u11  1)/At.  (A2) 

Also  define, 


ip  = du/dx  = u* . 


(A3) 


Then  Eq.  (Al)  can  be  written  as  the  equivalent  first-order  system 
given  by 

u'  = ip  , (A4) 

r = W + (u  - u”"1)  . (A5) 

If  Eqs.  (A4)  and  (A5)  are  viewed  as  characteristic  equations  for  an 
initial  value  problem,  then  characteristic  theory  states  that  the 
solution  of  Eqs.  (A4)  and  (AS)  generates  the  integral  surface  u(x,i(i)  of 
the  equation 

ux  + UiJ>  ^ + ~ un"1^  = + * (A6) 

The  general  solution  to  Eq.  (A6)  has  been  shown  by  Meyer  (Ref.  18)  to 
be 

u(x,iji)  = V(x)ip  + W(x)  . (A7) 


Meyer , G.  H. , Initial  Value  Methods  for  Boundary  Value  Problems 
(Theory  and  Application  of  Invariant  Imbedding) , 100th  Volume  in 
the  Series  Mathematics  in  Science  and  Engineering , ed.  R.  Bellman, , 
Academic  Press 3 NY  (1972). 


■45 


Substituting  Eq.  (A7)  into  Eq.  (A6)  yields  an  equation  for  both  V 


and  W, 

V'  = 1 - RV  - V2/At  , (A8) 

W'  = -V(W  - un_1)/At  . (A9) 

The  initial  conditions  required  for  these  two  equations  follow  from  a 
compliance  with  the  boundary  condition  (Alb)  at  -xq, 

V(-xq)  = 0,  W(-xo)  - 0 . , (A10) 

Substitution  of  Eq.  (A7)  into  Eq.  (A5)  gives  the  equation  governing  iJj, 

M = (R  + V/At)ip  + (W  - un"1)/At  . (All) 

The  initial  value  of  follows  directly  from  the  gradient  boundary 
condition  (Ale)  at  x=0,  i.e. 

MO)  = G | (A12) 


| new  time  level. 

The  sequence  of  computation  is  then  the  following: 

1.  Integrate  Eqs.  (A8)  and  (A9)  from  -x  to  0 based  on  the 
initial  values  given  by  Eq.  (A10) . 

2.  Satisfy  Eq.  (A12)  at  the  boundary  x=0. 

3.  Since  V(x)  and  W(x)  are  now  known,  integrate  Eq.  (All)  for 
ip(x)  from  0 to  -xq. 

4.  Evaluate  u(x)  at  the  new  time  level  from  Eq.  (A7) . 

It  may  be  noted  that  Eq.  (A8)  is  a Riccati  equation  which,  along 
with  the  initial  condition  Eq.  (A10) , has  the  exact  solution, 

V(x)  = [S  coth  e + R/2]"1  (A13) 

where  C=S(x+x)  [-x  < x < 0]  , 

o o 

and  S = [(R/2)2  + 1/At]1/2  . 

This  eliminates  the  need  to  numerically  integrate  Eq.  (A8)  in  the 
computation  sequence  listed  above. 
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Nonlinear  Equation 


Since  the  equation  and  the  gradient  boundary  condition  which 
govern  the  combustion  models  are  nonlinear,  a modification  to  the  above 
procedure  is  required.  For  these  models,  the  coefficient  R and  the 
function  G are  both  dependent  on  u(0)  at  the  new  time  level.  Hence  an 
iteration  cycle  is  created  between  steps  (1)  and  (2)  above;  (a)  an 
initial  guess  or  the  last  iterative  value  of  u(0)  is  used  to  integrate 
Eq.  (A8)  and  estimate  the  value  of  G,  (b)  Eq.  (A12)  then  determines 
^(0),  and  (c)  a new  value  of  u(0)  is  found  from  Eq.  (A7) . The  cycle 
must  be  repeated  until  convergence  is  Obtained.  The  computation  then 
continues  to  step  (3) . 
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